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Abstract 

We present an algorithm for the efficient samphng of conditional 
paths of stochastic differential equations (SDEs). While unconditional 
path sampling of SDEs is straightforward, albeit expensive for high di- 
mensional systems of SDEs, conditional path sampling can be difficult 
even for low dimensional systems. This is because we need to pro- 
duce sample paths of the SDE which respect both the dynamics of the 
SDE and the initial and endpoint conditions. The dynamics of a SDE 
are governed by the deterministic term (drift) and the stochastic term 
(noise) . Instead of producing conditional paths directly from the orig- 
inal SDE, one can consider a sequence of SDEs with modified drifts. 
The modified drifts should be chosen so that it is easier to produce sam- 
ple paths which satisfy the initial and endpoint conditions. Also, the 
sequence of modified drifts converges to the drift of the original SDE. 
We construct a simple Markov Chain Monte Carlo (MCMC) algorithm 
which samples, in sequence, conditional paths from the modified SDEs, 
by taking the last sampled path at each level of the sequence as an ini- 
tial condition for the sampling at the next level in the sequence. The 
algorithm can be thought of as a stochastic analog of deterministic ho- 
motopy methods for solving nonlinear algebraic equations or as a SDE 
generalization of simulated annealing. The algorithm is particularly 
suited for filtering/smoothing applications. We show how it can be 
used to improve the performance of particle filters. Numerical results 
for filtering of a stochastic differential equation are included. 



Introduction 

The study of systems arising in different areas, from signal processing and 
chemical kinetics to econometrics and finance (see e.g. ll^llSj) often requires 



1 



the sampling of paths of stochastic differential equations (SDEs) subject to 
initial and endpoint conditions. While unconditional path sampling of SDEs 
is straightforward, albeit expensive for high dimensional systems of SDEs, 
conditional path sampling can be difficult even for low dimensional systems. 
This is because we need to produce sample paths of the SDE which respect 
both the dynamics of the SDE and the initial and endpoint conditions. An 
analogous situation arises in ordinary differential equations, where it can be 
considerably more difficult to create solutions to boundary value problems 
than it is to construct solutions to initial value problems (see e.g. Ch. 8 in 
[1]). The problem of conditional path sampling of SDEs has been a subject 
of active research in recent years and some very interesting approaches have 
already been developed (see e.g. [2l ITSl [2T] ). 

The dynamics of a SDE are governed by the deterministic term (drift) 
and the stochastic term (noise). Instead of producing conditional paths 
directly from the original SDE, one can consider a sequence of SDEs with 
modified drifts. The modified drifts should be chosen so that it is easier 
to produce sample paths which satisfy the initial and endpoint conditions. 
Also, the sequence of modified drifts converges to the drift of the original 
SDE. We construct a simple Markov Chain Monte Carlo (MCMC) algorithm 
which samples, in sequence, conditional paths from the modified SDEs, by 
taking the last sampled path at each level of the sequence as an initial 
condition for the sampling at the next level in the sequence. 

We have used the drift relaxation algorithm to modify a popular filtering 
method called particle filter f5]. A particle filter is a sequential importance 
sampling algorithm which is based on the recursive (online) Bayesian up- 
dating of the values of samples (called particles) to incorporate information 
from noisy observations of the state of a dynamic model. While the par- 
ticle filter is a very versatile method it may require a very large number 
of samples to approximate accurately the conditional density of the state 
of the model. This has led to considerable research (see e.g. [3122]) into 
how one can modify a particle filter to make it more efficient (see also f3] 
for a different approach to particle filtering). As an application of the drift 
relaxation algorithm we show in Section [2] how it can be used to construct 
a more efficient particle filter. 

The paper is organized as follows. Section [T] presents the drift relaxation 
algorithm for an SDE conditional path sampling problem. Section [2] shows 
how to use the algorithm to modify a particle filter. Section [3] contains 
numerical results for the application of the modified particle filter to the 
standard example of filtering a diffusion in a double-well potential (more 
elaborate examples will be presented in [13] ) . Finally, Section [H discusses 
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the results as well as current and future work. 



1 Conditional path sampling and drift relaxation 

Suppose that we are given a system of stochastic differential equations 
(SDEs) 

dXt = a{Xt)dt + a{Xt)dBt, (1) 

Suppose also that we want to construct, in the time interval [0,T], sample 
paths from ([1]) such that the endpoints are distributed according to the 
densities h{Xo) and g{Xx) respectively. Equation ([1]) can be discretized in 
the interval [0, T] by some numerical approximation scheme [0]. Suppose 
that we have discretized the interval [0, T] using a stepsize At = T/I. To 
construct conditional paths of ([T]) we can start by sampling initial conditions 
from the density h. For each initial condition Xq we can then produce the 
desired conditional paths by sampling the density 

I 

1[p{Xt^Xt,_MXt), (2) 

i=l 

where p{X'j'.\Xx^_-i) is the transition probability from -^Ti_i at time Tj_i 
to the point Xt^ at time Tj. The density ([2]) can be sampled using MCMC 
assuming that the transition densities p(X7^ can be evaluated. How- 

ever, the major issue with the MCMC sampling is whether it can be per- 
formed efficiently (see e.g. [211 E]). Instead of MCMC sampling directly from 
the density ([2]) i.e., starting from an arbitrary initial path and modifying 
it to become a path corresponding to (l2|), we can aid the MCMC sampling 
process by providing the MCMC sampler of the density ([2]) with a better 
initial condition. 

To this end, consider an SDE system with modified drift 

dYt = b{Yt)dt + a{Yt)dBt, (3) 

where hiYt) can be suitably chosen to facilitate the conditional path sampling 
problem. 

Also, consider the collection of L -|- 1 modified SDE systems 

dY} = (1 - ei)b{Y})dt + €ia{Y})dt + (T{Yl)dBt, 

where G [0, 1], / = 0, . . . , L, with e/ < e^+i, eo = and = 1. Instead of 
sampling directly a (conditional) path from the SDE ([1]), one can sample a 
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path from the modified SDE ([3]) and graduaUy morph the path into a path 
of dH). 

Drift relaxation algorithm: 

• (/ = 0) Begin with a sample path from the modified SDE ([3]). 

• Sample through MCMC the density nLiP(^Tj^T,_i)5(^r)- 

• For / = 1, ...,L take the last sample path from the (/ — l)st SDE and 
use it as in initial condition for MCMC sampling the density 

/ 

i=l 

at the Ith level. 

• Keep the last sample path at the Lth level. 

The drift relaxation algorithm is similar to Simulated Annealing (SA) used 
in equilibrium statistical mechanics [TO]. However, instead of modifying a 
temperature as in SA, here we modify the drift of the system. 

2 Application to particle filtering 

We show in this section how the drift relaxation algorithm can be applied 
to particle filtering with the aim of bringing the samples closer to the ob- 
servations. 

2.1 Generic particle filter 

Suppose that we are given an SDE system and that we also have access 
to noisy observations Zti , ■ ■ ■ , Zt/^ of the state of the system at specified 
instants Ti, . . . ,Tk- The observations are functions of the state of the sys- 
tem, say given by Zj-^. = G{XTf.,£,k), where ^fc,A; = 1,...,K are mutually 
independent random variables. For simplicity, let us assume that the distri- 
bution of the observations admits a density g{XTf., Zt^), i.e., p{Zt^\Xt^) oc 
g{XT„ZT,). 

The filtering problem consists of computing estimates of the conditional 
expectation E[f{XT^)\{ZTj}j^i\, i.e., the conditional expectation of the 
state of the system given the (noisy) observations. Equivalently, we are 
looking to compute the conditional density of the state of the system given 
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the observations p{XT^.\{ZTj}j=i)- There are several ways to compute this 
conditional density and the associated conditional expectation but for prac- 
tical applications they are rather expensive. 

Particle filters fall in the category of importance sampling methods. Be- 
cause computing averages with respect to the conditional density involves 
the sampling of the conditional density which can be difficult, importance 
sampling methods proceed by sampling a reference density q{XTf,\{ZTj}j=i) 
which can be easily sampled and then compute the weighted sample mean 

^ 1 J!^, p{X^ |{^T }j^=i) 



or the related estimate 



q{X-\{ZTM^) 



where has been replaced by the approximation 

Particle filtering is a recursive implementation of the importance sampling 
approach. It is based on the recursion 

p{Xt,\{Zt^}]=^) a g{XT,,ZTMXTMZT,]]zl), (5) 
where p(XTj{ZT,}y) = |p(XtJXt,_ Jp(Xt,_ J{Zt, )dXT,_, . (6) 

If we set 

q{XTMZT^}U)=piXn\{ZT,})zl), 

then from ([5]) we get 

p{XtMZt^}U) ,^ 7 , 
.fc N 5(^rfc , Zt^ ) . 

g(^Tj{^T,}j=l) 

The approximation in expression @ becomes 

)5(^r ,^tJ 

l^n=l9\Xj' .ZtJ 
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From ([7|) we see that if we can construct samples from the predictive distri- 
bution p{XT^.\{ZTj }j=i ) then we can define the (normalized) weights = 

-, use them to weigh the samples and the weighted samples 



will be distributed according to the posterior distribution p(^Tfc |{-^Tj }j=i)- 
In many applications, most samples will have a negligible weight with 
respect to the observation, so carrying them along does not contribute sig- 
nificantly to the conditional expectation estimate (this is the problem of 
degeneracy |10|). To create larger diversity one can resample the weights to 
create more copies of the samples with significant weights. The particle filter 
with resampling is summarized in the following algorithm due to Gordon et 
al. [S]. 

Particle filter 

1. Begin with unweighted samples X'^^ ^ from p{Xti._j^\{Ztj}^zI) ■ 

2. Prediction: Generate samples from p(Xrj, 

3. Update: Evaluate the weights 

g(Xlp , Zt, ) 
' En=i 9{X!P.ZtJ 



4. Resampling: Generate N independent uniform random variables 
{0'^}^^i in (0, 1). Forn = 1, . . . , TV let XJ^ = where 

1=1 1=1 
where j can range from 1 to N. 

5. Set k = k + \ and proceed to Step 1. 



The particle filter algorithm is easy to implement and adapt for dif- 
ferent problems since the only part of the algorithm that depends on the 
specific dynamics of the problem is the prediction step. This has led to the 
particle filter algorithm's increased popularity [5]. However, even with the 
resampling step, the particle filter can still need a lot of samples in order to 
describe accurately the conditional density p[XT,.\{ZTj}^=i) ■ Snyder et al. 
[I6j have shown how the particle filter can fail in simple high dimensional 
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problems because one sample dominates the weight distribution. The rest 
of the samples are not in statistically significant regions. Even worse, as we 
will show in the numerical results section, there are simple examples where 
not even one sample is in a statistically significant region. In the next sub- 
section we present how drift relaxation can be used to push samples closer 
to statistically significant regions. 

2.2 Particle filter with MCMC step 

Several authors (see e.g. |7i [22]) have suggested the use of a MCMC step 
after the resampling step (Step 4) in order to move samples away from 
statistically insignificant regions. There are many possible ways to append 
an MCMC step after the resampling step in order to achieve that objective. 
The important point is that the MCMC step must preserve the conditional 
density p(^rfc K-^T, }j=i)- In the current section we show that the MCMC 
step constitutes a case of conditional path sampling. 

We begin by noting that one can use the resampling step (Step 4) in the 
particle filter algorithm to create more copies not only of the good samples 
according to the observation, but also of the values (initial conditions) of the 
samples at the previous observation. These values are the ones who have 
evolved into good samples for the current observation (see more details in 
[22]). The motivation behind producing more copies of the pairs of initial 
and final conditions is to use the good initial conditions as starting points 
to produce statistically more significant samples according to the current 
observation. This process can be accomplished in two steps. First, Step 4 
of the particle filter algorithm is replaced by 

Resampling: Generate independent uniform random variables {0""}^ 
in (0, 1). Forn = 1, . . . ,7V let (X^^_^,XJJ = (X^^_^,X^Jwhere 

1=1 1=1 

Also, through Bayes' rule [22] one can show that the posterior density 
p{XT,^\{ZTj}^j^i) is preserved if one samples from the density 

^(^Tfe , ^Tjp(^Tj^Tfc_i ) 

where Xt^.^^ are given by the modified resampling step. This is a problem 
of conditional path sampling for (continuous-time or discrete) stochastic 
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systems. The important issue is to perform the necessary sampHng efficiently 
[21 [22]. We propose to do that here using drift relaxation (see Section [TJ. 
The particle filter with MCMC step algorithm is given by 

Particle filter with MCMC step 

1. Begin with unweighted samples XJj^^ ^ from p{XTf,_^\{ZTj}^jZ\) ■ 

2. Prediction: Generate samples X'^^ ivom. p{Xt^\Xt^_^). 

3. Update: Evaluate the weights 

4. Resampling: Generate N independent uniform random variables 

in (0,1). For n = 1,...,N let (X^,_,,X- ) = (X^,_^,X^J 

where 

1=1 1=1 

where j can range from 1 to N. 

5. MCMC step: For n = 1,...,N choose a modified drift (possibly 
different for each n). Construct a path for the SDE with the modified 
drift starting from X^^ ^ . Construct through drift relaxation a Markov 
chain for with stationary distribution 

6. Set X" = . 

7. Set k = k + 1 and proceed to Step 1. 

Since the samples X^^ = 1^ are constructed by starting from different 
sample paths, they are independent. Also, note that the samples X^^ are 
unweighted. However, we can still measure how well these samples approx- 
imate the posterior density by comparing the effective sample sizes of the 
particle filter with and without the MCMC step. For a collection of A^ 
samples the effective sample size ess{Tk) is defined by 
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where 



^E(5(^T,>^Tj-^fc)2 and 1^, = 1^5(X^^,ZtJ. 

n=l n=l 

The effective sample size can be interpreted as that the N weighted samples 
are worth of ess{Tk) = _/^2 i-i-d. samples drawn from the target density, 

which in our case is the posterior density. By definition, ess{Tk) < N. If the 
samples have uniform weights, then ess{Tk) = N. On the other hand, if all 
samples but one have zero weights, then ess{Tk) = 1. 



3 Numerical results 

We present numerical results of the particle filter algorithm with MCMC 
step for the standard problem of diff'usion in a double-well potential (more 
elaborate applications of the method will be presented elsewhere [l3]). Our 
objective here is to show how the generic particle filter's performance can be 
significantly improved by incorporating the MCMC step via drift relaxation. 

The problem of diffusion in a double well potential is described by the 
scalar SDE 

dXt = -4.Xt{Xf - 1) + ^Bt. (8) 

The deterministic part (drift) describes a gradient fiow for the potential 
U (x) = x^ — which has two minima, at x = ±1. In the notation of Section 
[2] we have a{Xt) = —4Xt{Xf — 1) and (j{Xt) = ^. If the stochastic term is 
zero the solution wanders around one of the minima depending on the value 
of the initial condition. A weak stochastic term leads to rare transitions 
between the minima of the potential. We have chosen the coefficient ^ to 
make the stochastic term rather weak. This is done because we plan to 
enforce the observations to alternate among the minima, and thus check if 
the particle filter can track these transitions. 

The SDE ^ is discretized by the Euler-Maruyama [9] scheme with step 
size At = 10~^ which is small enough to guarantee stability of the scheme. 
The initial condition is set to —1 and there is a total of 10 observations at 
Tk = k,k = 1, ... ,10. The observations are given by Zxf. = Xt^ + ^fc, where 

~ N(0, 0.01) for k = 1, . . . , 10. For this choice of observation noise, the 
observation density (also called likelihood) is given by 



9{Xt^,ZtJ oc exp 



2*0.01 



(9) 
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The observations alternate between 1 and —1. In particular, for A; = 1, . . . , 10 
we have Zt^. = — 1 if /c is odd and Zt^. = 1 if A; is even. 

In order to apply the MCMC step with drift relaxation we need to define 
the modified drift h{Yt) for the process Yt given by 

dYt = b{Yt) + ^dBt. (10) 

The modified drift can be the same for all the samples or different for each 
sample. Since the difficulty in tracking the observations comes from the 
inability of the original SDE ([8]) to make frequent transitions between the 
two minima of the double well, an intuitively appealing choice for b(Yf) 
is b{Yt) = —aAYt{Y^ — 1), where a < 1. This drift corresponds to the 
potential W{y) = a{y'^ — 2y^). The potential W{y) has its minima also 
located at y = ±1. However, the value of the potential at the minima is —a 
instead of —1 for the potential U{x). This means that the wells corresponding 
to the minima of W{y) are shallower than the wells corresponding to the 
minima of U{x). This makes the transitions between the two wells for the 
process Yt more frequent than for the original process Xt- For the numerical 
experiments we have chosen a = 0.1. 

The sequence of modified SDEs for the drift relaxation algorithm with 
L levels is given by 

dYl = (1 - ei)b{Yl)dt + eia{Yl)dt + ]^dBu (11) 

where e; G [0, 1], / = 0, . . . ,-L, with ei < e^+i, eo = and €l = 1- For our 
numerical experiments we chose L = 10 and e/ = I/IO. 

Recall that the density we want to sample during the MCMC step is 
given by 

g{XT„ZTMXn\XT,_,), 

where p{Xt^\Xt^_-^) is the transition probability between Xti._^ and ^r^- 
For many applications, sampling directly from piX^f, \Xti^_-^) may be impos- 
sible. Thus, one needs to resort to some numerical approximation scheme 
which approximates the path between Xt,._-^ and Xt^. by a discretized path. 
However (see [22] for details), even the evaluation of the discretized path's 
density may not be efficient. Instead, by using the fact that each Brown- 
ian path in ([8j) gives rise to a unique path for Xt [E], we can replace the 
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sampling of g{XT^, ZTf.)p{XT^\XT^_^) by sampling from the density 



exp 



(Z, 



2*0.01 



exp 



(AS, 



n^2^ 

i 



2* At 



2 * 0.01 



7-1 



2* At 



n\2 \ 1 



(12) 



where {ABf}^^Q are the Brownian increments of the discretized path con- 
necting XTf._-^ and Xxf.- Also, note that the final point Xt,. has now become 
a function of the entire Brownian path {ABf}l~Q . For the numerical exper- 
iments we have chosen At = = 10~ which, since — T^-i = 1, 

gives / = 100. 

We use drift relaxation to produce samples from the density (fT2]l . The 
Markov chain at each level of the drift relaxation algorithm is constructed 
using Hybrid Monte Carlo (HMC) [10]. At the Ith level, we can discretize 
pip, say with the Euler-Maruyama scheme, and the points on the path will 
be given by 



^/)K>^(Si)A,)Ai + e,aiYl:-_,^^^)At + Ub 



Ln 



^1' 



for i = 1,...,/. We can use more sophisticated schemes than the Euler- 
Maruyama scheme for the discretization of the simplified SDE (jlOp at the 
cost of making the expression for the density more complicated. 

We can define a potential ({Ai?^' }^_q) for the variables {Ai?j' }^_q . 
The potential is given by 



K,({{AS, }.=o)- j:^^ + 2^ 



i=0 . ^ ^ _ 

' 2* At 



i=0 



and the density to be sampled can be written as 



exp 



The subscript ei is to denote the dependence of the potential on the drift 
relaxation parameter e;. In HMC one considers the variables on which the 
potential depends as the position variables of a Hamiltonian system. In our 
case we have / position variables so we can define a /-dimensional position 
vector {qiYi^i- The next step is to augment the position variables vector by 
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a vector of associated momenta {pi}i=i- Together they form a Hamiltonian 
system with Hamiltonian given by 



T 
p p 



where p = {pi, . . . ,pi) is the momenta vector. Thus, the momenta variables 
are Gaussian distributed random variables with mean zero and variance 1. 
The equations of motion for this Hamiltonian system are given by Hamilton's 
equations 

^ = ^ and ^ = -^ for ^ = 1,...,I. 

dr opi dr oqi 

Note that the Hamiltonian depends also on the initial condition for each 
sample Yq' and we have written an equation for the evolution of qi = Yq" 
as well its associated momentum pi. Since the Y^" are fixed by the resam- 
pling procedure we do not evolve them. However, note that the Brownian 
increment ABq needs to be evolved because it affects the evolution of Ya*- 
HMC proceeds by assigning initial conditions to the momenta variables 

T 

(through sampling from exp( — 2^)), evolving the Hamiltonian system in 
fictitious time r for a given number of steps of size 6t and then using the 
solution of the system to perform a Metropolis accept/reject step (more de- 
tails in pO]). After the Metropolis step, the momenta values are discarded. 
The most popular method for solving the Hamiltonian system, which is the 
one we also used, is the Verlet leapfrog scheme. In our numerical imple- 
mentation, we did not attempt to optimize the performance of the HMC 
algorithm. For the sampling at each level of the drift relaxation process we 
used 10 Metropolis accept/reject steps and 1 HMC step of size 6t = 10~^ 
to construct a trial path. A detailed study of the drift relaxation/HMC 
algorithm for conditional path sampling problems outside of the context of 
particle filtering will be presented in a future publication. 

For the chosen values of the parameters for the drift relaxation and 
HMC steps, the particle filter with MCMC step is about 500 times more 
expensive per sample (particle) than the generic particle filter. However, we 
show that this increase in cost per sample is worthwhile. Figure [J compares 
the performance of the particle filter with MCMC step with 10 samples and 
the generic particle filter with 5000 samples. It is obvious that the particle 
filter with MCMC step follows accurately all the transitions between the 
two minima of the double-well. On the other hand, the generic particle 
filter captures accurately only every other observation. It fails to performs 
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Figure 1: Comparison of the conditional expectation of Xt as computed by 
the generic particle filter and the particle filter with MCMC step. Different 
drift for each sample. 



the transitions between the two minima of the double-well. Of course, since 
the particle filter with MCMC step uses only 10 samples the conditional 
expectation estimate of the hidden signal is not as smooth as the estimate 
of the generic particle filter which uses 5000 samples. However, the particle 
filter with MCMC step allows a much better resolution of the conditional 
density (conditioned on the observations). This can be seen by computing 
the effective sample size for the two filters. Figure [2] shows the effective 
sample size for both filters. Because of the different number of samples used 
in the two filters we have plotted the effective sample size for each filter as a 
percentage of the respective sample size. We see that the particle filter with 
MCMC step has overall a much better sample size than the generic particle 
filter. The generic particle filter has a wildly fiuctuating effective sample size. 
In particular, since the generic particle filter misses every other observation, 
the corresponding effective sample size dips down to 1 sample for the missed 
observations. Note from Figured] that even this one sample which dominates 
the observation weight does not come close to the observation. For the 
observations that the generic particle filter does capture, its effective sample 
size is still lower than the effective sample size of the particle filter with 
MCMC step. 
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«-# Effective sample size for PF with MCMC (% of total samples) 
#—9 Effeclive sample size for PF {% of total samples) 




Observations 



Figure 2: Comparison of the effective sample size as computed by the generic 
particle filter and the particle filter with MCMC step. A different drift is 
used for each sample. 

4 Discussion 

We have presented an algorithm for conditional path sampling of SDEs. 
The proposed algorithm is based on drift relaxation which allows to sample 
conditional paths from a modified drift equation. The conditional paths of 
the modified drift equation are then morphed into conditional paths of the 
original equation. We have called this process of gradually enforcing the drift 
of the original equation drift relaxation. The algorithm has been used to 
create a modified particle filter for SDEs. We have shown that the modified 
particle filter's performance is significantly better than the performance of 
a generic particle filter. 

In the current work, we have examined the application of drift relax- 
ation to the filtering problem of diffusion in a double-well potential which 
is a standard example in the filtering literature. The same algorithm can 
be applied to the problem of tracking a single target. A problem of great 
practical interest is that of tracking not only one but multiple moving tar- 
gets [T^ [T9\ [20] . The multi-target tracking problem is much more dif- 
ficult than the single-target problem due to the combinatorial explosion of 
the number of possible target-observation association arrangements. In this 
context, the accurate tracking of each target becomes crucial. Suppose that 
only one of the targets is of interest and the rest act as decoys |12j . The in- 
ability to track each potential target accurately can lead to ambiguity about 
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the targets' movement if the observations for different targets are close. We 
have already applied the drift relaxation modified particle filter to multi- 
target tracking problems with very encouraging results which will appear 
elsewhere [13]. 
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